1 Introduction

Contexte :

Nous somme employé d’Enercoop, jeune entreprise spécialisé dans la fourniture d’énergie renouvelable

But de l’étude :

  • Prédire la consommation électrique sur 1 an

  • Crée un modèle de prévision

Pour cela, nous allons utilisé plusieurs modèles, plusieurs méthodes, de prédiction notamment :

  • Lissage exponentielle

  • SARIMA

  • Machine learning

  • Prophet

Bibliographie :

## Install library
library(tidyverse)
library(lubridate)
library(fabletools)
library(fable)
library(tsibble)
library(feasts)
library(car)
library(dygraphs)
library(broom)

2 Imports des données

Nous restreignons l’analyse à la région Grand Est avec comme capteur Strasbourg Entzheim.
La période s’étend 2013 et 2019

2.1 Consommation

Source https://www.rte-france.com/eco2mix/telecharger-les-indicateurs

2.1.1 Imports des consommations electrique

2.1.2 Création du df

## Ajout de toutes les années
data_conso_grandest <- rbind(RTE2019,RTE2018,RTE2017,RTE2016,RTE2015,RTE2014,RTE2013)
## modification du df
conso_daily_grandest <-  data_conso_grandest %>% 
                        filter(Consommation > 0) %>% # suppression des valeurs a 0
                        mutate(Date_only = ymd(Date)) %>% # conversion en year month day
                        mutate(Date_only = as_date(Date_only)) %>% # conversion au format date
                        mutate(Renouvelable = (Eolien+Hydraulique+Solaire+Pompage+Bioénergies))  %>% # Creation du variable energie renouvelable
                        select(Date_only, Consommation,Renouvelable) %>% # Selection des variables d'interet
                        group_by(Date_only) %>% # group by daily avec moyenne 
                        summarise_all(funs(mean(.,na.rm=TRUE))) %>% 
                        as_tsibble(index= Date_only) %>%  # conversion du df en tsibble
                        ungroup()
## Creation d'un dataframe month_avg
conso_month_grandest <- conso_daily_grandest %>% 
                            index_by(yearmonth(Date_only)) %>% # group_by(index_by) month
                            summarise_all(funs(mean(.,na.rm=TRUE)))

2.1.3 Analyse du df

DataExplorer::plot_intro(conso_daily_grandest)

questionr::describe(conso_daily_grandest)
## [2556 obs. x 3 variables] tbl_ts tbl_df tbl data.frame
## 
## $Date_only: 
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
## 
## $Consommation: 
## numeric: 4672.10416666667 5591.9375 6086.375 5956.52083333333 5408.29166666667 5086.875 6042.875 6619.3125 6852.70833333333 6856.6875 ...
## min: 3295.33333333333 - max: 7849.14583333333 - NAs: 0 (0%) - 2534 unique values
## 
## $Renouvelable: 
## numeric: 2083.75 1396.75 1644.20833333333 1317.85416666667 1225.04166666667 1166.625 1166.29166666667 1133.875 1163.5625 1171.91666666667 ...
## min: 349.270833333333 - max: 4020.5625 - NAs: 0 (0%) - 2522 unique values
conso_daily_grandest %>% 
  psych::pairs.panels(
              method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE # show correlation ellipses
             )

Pas de données manquante ou de ligne incomplète
2556 observations = 7 ans
La consommation au fil du temps a une moyenne stable = saisonnalité additive

2.2 Température

2.2.1 Imports des températures relevés

Source : https://www.ncdc.noaa.gov/cdo-web/search

2.2.2 Création du df

temp_strasbourg_2013_2019 <- mutate(temp_strasbourg_2013_2019,DATE = ymd(DATE)) %>% # conversion en year month day
                            select(-SNWD) %>% # suppression de la colonne snow
                            as_tsibble(index="DATE") # conversion du df en tsibble

2.2.3 Analyse du df

temp_strasbourg_2013_2019 %>% questionr::describe() 
## [2556 obs. x 6 variables] tbl_ts tbl_df tbl data.frame
## 
## $NAME: 
## character: "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" ...
## NAs: 0 (0%) - 1 unique values
## 
## $DATE: 
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
## 
## $PRCP: 
## numeric: 1 0 0 0.8 0 0 0 0 0.4 0.6 ...
## min: 0 - max: 66.3 - NAs: 0 (0%) - 153 unique values
## 
## $TAVG: 
## numeric: 6.6 3.6 3.2 7.1 7.2 6.7 4.7 3.4 2.2 2.9 ...
## min: -7.9 - max: 29.6 - NAs: 0 (0%) - 324 unique values
## 
## $TMAX: 
## numeric: 8.9 7.2 6.1 8.7 7.9 8.4 9.2 3.8 2.5 4.2 ...
## min: -4.9 - max: 38.9 - NAs: 0 (0%) - 396 unique values
## 
## $TMIN: 
## numeric: 3.8 0.3 -0.8 5.2 6.3 6 0.2 2.2 1.9 1.3 ...
## min: -11.6 - max: 23.8 - NAs: 0 (0%) - 290 unique values
temp_strasbourg_2013_2019 %>% DataExplorer::plot_intro()

temp_strasbourg_2013_2019 %>% psych::pairs.panels(
                                          method = "pearson", # correlation method
                                         hist.col = "#00AFBB",
                                         density = TRUE,  # show density plots
                                         ellipses = TRUE # show correlation ellipses
                                         )

Ici aussi, pas de problème particulier dans les données

2.3 DJU ( degrés jour unifiés)

2.3.1 Imports des DJU

Source : https://cegibat.grdf.fr/simulateur/calcul-dju

## Imports données (Grand est 67)
dju_month_grandest <- read.csv("data/calcul_DJU_05_05_2021.csv",dec=",")

2.3.2 Modification du df

dju_month_grandest <- dju_month_grandest %>% 
                      mutate (date = ym(DATE)) %>% 
                      relocate(DJU, .after = last_col()) %>% 
                      as_tsibble(index = date)

2.3.3 Création d’un data daily par calcul

En utilisation la méthode de calcul Méteo de Cegibat

seuil_temp = 18 # seuil point de changement température
dju_daily_grandest <- temp_strasbourg_2013_2019 %>% 
                      mutate(
                      dju_calc_chauffage = case_when( ## calcul DJU chauffage
                      seuil_temp <= ((TMIN + TMAX)/2) ~ 0,
                      seuil_temp > ((TMIN + TMAX)/2) ~ seuil_temp - ((TMIN + TMAX)/2),
                      TRUE ~ 0),
                      dju_calc_clim = case_when( ## Calcul DJU climatisation
                      seuil_temp >= ((TMIN + TMAX)/2) ~ 0,
                      seuil_temp < ((TMIN + TMAX)/2) ~ ((TMIN + TMAX)/2) -seuil_temp,
                      TRUE ~ 0 ),
                      dju_calc_total = dju_calc_chauffage + dju_calc_clim # calcul DJU total
                      )

2.3.4 Création df monthly

dju_month_grandest_calc <- dju_daily_grandest %>% 
                            index_by(yearmonth(DATE)) %>% 
                            summarise_at(vars(dju_calc_chauffage:dju_calc_total),~sum(.))

2.3.5 Comparaision calcul vs cegibat

autoplot(dju_month_grandest_calc,dju_calc_chauffage,color='gray')+
  geom_line(aes(y=dju_month_grandest$DJU), colour = "#D55E00")+
  labs(
    y = "DJU",
    title = "Comparaision calcul vs cegibat"
  )

Les graphiques sont proches, nous pouvons valider la méthode de calcul

3 Exploration des données

3.1 Création des df finaux

## merge des dataframmes
data_est_day <- left_join(conso_daily_grandest, dju_daily_grandest, by = c('Date_only'='DATE')) %>% 
  select(-NAME) %>%  # Suppresion variable name 
  as_tsibble(index = Date_only)
## Creation d'un data month
data_est_month <- data_est_day %>% # group_by
  index_by(yearmonth(Date_only)) %>% 
  summarise(
  across(c(TAVG:TMIN),~mean(.)),
  across(c(Consommation:Renouvelable,dju_calc_chauffage:dju_calc_total,PRCP),~sum(.))
            )

3.2 Analyse des DF

data_est_day %>% questionr::describe() 
## [2556 obs. x 10 variables] tbl_ts tbl_df tbl data.frame
## 
## $Date_only: 
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
## 
## $Consommation: 
## numeric: 4672.10416666667 5591.9375 6086.375 5956.52083333333 5408.29166666667 5086.875 6042.875 6619.3125 6852.70833333333 6856.6875 ...
## min: 3295.33333333333 - max: 7849.14583333333 - NAs: 0 (0%) - 2534 unique values
## 
## $Renouvelable: 
## numeric: 2083.75 1396.75 1644.20833333333 1317.85416666667 1225.04166666667 1166.625 1166.29166666667 1133.875 1163.5625 1171.91666666667 ...
## min: 349.270833333333 - max: 4020.5625 - NAs: 0 (0%) - 2522 unique values
## 
## $PRCP: 
## numeric: 1 0 0 0.8 0 0 0 0 0.4 0.6 ...
## min: 0 - max: 66.3 - NAs: 0 (0%) - 153 unique values
## 
## $TAVG: 
## numeric: 6.6 3.6 3.2 7.1 7.2 6.7 4.7 3.4 2.2 2.9 ...
## min: -7.9 - max: 29.6 - NAs: 0 (0%) - 324 unique values
## 
## $TMAX: 
## numeric: 8.9 7.2 6.1 8.7 7.9 8.4 9.2 3.8 2.5 4.2 ...
## min: -4.9 - max: 38.9 - NAs: 0 (0%) - 396 unique values
## 
## $TMIN: 
## numeric: 3.8 0.3 -0.8 5.2 6.3 6 0.2 2.2 1.9 1.3 ...
## min: -11.6 - max: 23.8 - NAs: 0 (0%) - 290 unique values
## 
## $dju_calc_chauffage: 
## numeric: 11.65 14.25 15.35 11.05 10.9 10.8 13.3 15 15.8 15.25 ...
## min: 0 - max: 26.25 - NAs: 0 (0%) - 548 unique values
## 
## $dju_calc_clim: 
## numeric: 0 0 0 0 0 0 0 0 0 0 ...
## min: 0 - max: 11.2 - NAs: 0 (0%) - 226 unique values
## 
## $dju_calc_total: 
## numeric: 11.65 14.25 15.35 11.05 10.9 10.8 13.3 15 15.8 15.25 ...
## min: 0 - max: 26.25 - NAs: 0 (0%) - 603 unique values
data_est_day %>% DataExplorer::plot_intro()

data_est_day %>% psych::pairs.panels(
                                          method = "pearson", # correlation method
                                         hist.col = "#00AFBB",
                                         density = TRUE,  # show density plots
                                         ellipses = TRUE # show correlation ellipses
                                         )

3.2.1 Analyse graphique de la variable d’intéret consommation

autoplot(data_est_day,Consommation)

La courbe ci dessus présente une saisonnalité annuelle avec une tendance stable

data_est_day %>% gg_season(Consommation, period = "year", labels = "both")  +
  labs(y="MW", title="Electricity demand: Grand Est")

data_est_day %>% gg_season(Consommation, period = "week", labels = "both")  +
  labs(y="MW", title="Electricity demand: Grand Est")

ggplot(data_est_day, aes(y=Consommation,x=TAVG))+
    geom_point()+
    theme_classic()+
    geom_smooth(method="lm", 
                colour="red", 
                formula=y~x+I(x^2))

Nous voyons bien ici une corrélation polynomiale de type y = x + x² entre la température et la consommation électrique

## Visualisation 
scatterplot(Consommation~dju_calc_total,data=data_est_day)

Idem pour ce graphique, nous notons une relation linéaire entre la consommation et les DJU

3.3 Correction de l’effet température

Ici, nous allons crée une variable consommation corrigé grâce aux régressions

3.3.1 Correction avec DJU

3.3.1.1 Régression linéaire

lm_temp_corr<-lm(data_est_day$Consommation~data_est_day$dju_calc_total)
summary(lm_temp_corr)
## 
## Call:
## lm(formula = data_est_day$Consommation ~ data_est_day$dju_calc_total)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1818.6  -406.6   102.1   427.4  1504.1 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4166.502     18.989  219.42   <2e-16 ***
## data_est_day$dju_calc_total  141.029      2.004   70.37   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 557.5 on 2554 degrees of freedom
## Multiple R-squared:  0.6597, Adjusted R-squared:  0.6596 
## F-statistic:  4952 on 1 and 2554 DF,  p-value: < 2.2e-16
a<-coef(lm_temp_corr)[2]
b<-coef(lm_temp_corr)[1]

3.3.1.2 Vérification des conditions de validité

Nous vérifions maintenant les conditions de validité de la régression linéaire

## Independance des résidus
acf(residuals(lm_temp_corr), main="lm_temp_corr")

durbinWatsonTest(lm_temp_corr)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6127372     0.7718504       0
##  Alternative hypothesis: rho != 0
# H0 (null hypothesis): There is no correlation among the residuals.
# HA (alternative hypothesis): The residuals are autocorrelated.

Nous voyons une forte autocorrelation des résidus sur le graphe et le test de Durbin Watson nous pousse à rejeter l’hypothèse de non corrélation

## Evaluation de l’hypothèse de normalité des résidus
plot(lm_temp_corr,2)

shapiro.test(residuals(lm_temp_corr))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm_temp_corr)
## W = 0.97373, p-value < 2.2e-16

La normalité des résidus est accepté

## Evaluation de l’hypothèse d’homogénéité des résidus
plot(lm_temp_corr, 3)

ncvTest(lm_temp_corr)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.916416, Df = 1, p = 0.047817

L’homoscédasticité des residus est accepté
Les conditions de validité du modèle sont validées

3.3.1.3 Consommation révisé

data_est_day <-
  data_est_day %>% mutate(conso_corrige_dju = Consommation - a * dju_calc_total)
data_est_month <-
  data_est_month %>% mutate(conso_corrige_dju = Consommation - a * dju_calc_total)
library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
data <- xts(x=data_est_day$conso_corrige_dju,order.by = data_est_day$Date_only)
dygraph(data)%>% dyRangeSelector()
data_est_month %>% 
autoplot(conso_corrige_dju)

3.3.2 Correction avec température

3.3.2.1 Régression polynomiale

lm_temp_corr2<-lm(data_est_day$Consommation~data_est_day$TAVG+I(data_est_day$TAVG^2))
summary(lm_temp_corr2)
## 
## Call:
## lm(formula = data_est_day$Consommation ~ data_est_day$TAVG + 
##     I(data_est_day$TAVG^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1787.4  -414.3   110.0   387.3  1575.6 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6793.0822    24.7485  274.49   <2e-16 ***
## data_est_day$TAVG      -189.4609     4.4057  -43.00   <2e-16 ***
## I(data_est_day$TAVG^2)    3.5738     0.1765   20.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 524.6 on 2553 degrees of freedom
## Multiple R-squared:  0.6989, Adjusted R-squared:  0.6986 
## F-statistic:  2963 on 2 and 2553 DF,  p-value: < 2.2e-16
x<-coef(lm_temp_corr2)[2]
x2<-coef(lm_temp_corr2)[3]

3.3.2.2 Vérification des conditions de validité

Nous vérifions maintenant les conditions de validité de la régression linéaire

## Indépendance  des residus
acf(residuals(lm_temp_corr2), main="lm_temp_corr")

durbinWatsonTest(lm_temp_corr2)
##  lag Autocorrelation D-W Statistic p-value
##    1        0.594772     0.8073433       0
##  Alternative hypothesis: rho != 0
# H0 (null hypothesis): There is no correlation among the residuals.
# HA (alternative hypothesis): The residuals are autocorrelated.

Nous voyons une forte autocorrelation des résidus sur le graphe et le test de Durbin Watson nous pousse à rejeter l’hypothèse de non corr&lation

## Evaluation de l’hypothèse de normalité des résidus
plot(lm_temp_corr2,2)

shapiro.test(residuals(lm_temp_corr2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm_temp_corr2)
## W = 0.97537, p-value < 2.2e-16

La normalité des résidus est accepté

## Evaluation de l’hypothèse d’homogénéité des résidus
plot(lm_temp_corr2, 3)

ncvTest(lm_temp_corr2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 16.43193, Df = 1, p = 5.0429e-05

L’homoscédasticité des résidus est accepté
Les conditions de validité du modèle sont validé

3.3.2.3 Consommation révisé

data_est_day <-
  data_est_day %>% mutate(conso_corrige_temp = Consommation - abs((x * TAVG + (x2^2)*TAVG)))
data_temp <- data_est_day %>% 
  index_by(yearmonth(Date_only)) %>% 
  summarise(sum(conso_corrige_temp))
## merge des dataframmes
data_est_month <- left_join(data_est_month, data_temp, by = 'yearmonth(Date_only)') 
library(xts)
data <- xts(x=data_est_day$conso_corrige_temp,order.by =data_est_day$Date_only)
dygraph(data)%>% dyRangeSelector()

3.4 Ajustements et transformation

Afin d’améliorer la compréhention et la précision du modèle, nous pouvons transformer celui ci.Soit en désaisonnalisant ou en appliquant une transformation de box-cox

3.4.1 Desaisonalisation Moving average

Il s’agit d’une des premières méthode de desaisonnalisation, utilisé depuis 1920. Pour estimer la tendance mensuel avec une saisonnalité annuel, nous utilisons une moyenne mobile 2x12MA

data_est_month_ma <- data_est_month %>%
  mutate(
    `12-MA` = slider::slide_dbl(Consommation, mean,
                .before = 5, .after = 6, .complete = TRUE),
    `2x12-MA` = slider::slide_dbl(`12-MA`, mean,
                .before = 1, .after = 0, .complete = TRUE)
  )
data_est_month_ma %>%
  autoplot(Consommation, colour = "gray") +
  geom_line(aes(y = `2x12-MA`), colour = "#D55E00") +
  labs(y = "MW",
       title = "Consommation")

3.4.2 Desaisonalisation STL

Une autre méthode utilisé par le census bureau est la méthode STL pour les séries additive sans variation calendaire, son principal avantage est de pouvoir traiter tout type de saisonnalité (heures, jours, mois, autres… )

3.4.2.1 Data month

## Desaisonnalisation de Consommation avec STL 
dcmp_stl_conso <- data_est_month %>%
  model(stl = STL(Consommation))
components(dcmp_stl_conso) %>%
  as_tsibble() %>%
  autoplot(Consommation, color="gray") +
  geom_line(aes(y=trend+remainder), color = "#D55E00") +
  labs(
    y = "MV",
    title = "electricity consommation"
  )

components(dcmp_stl_conso) %>% autoplot()

## Desaisonnalisation de TAVG
dcmp_stl_tavg <- data_est_month %>%
  model(stl = STL(TAVG))
components(dcmp_stl_tavg) %>%
  as_tsibble() %>%
  autoplot(TAVG, color="gray") +
  geom_line(aes(y=trend), color = "#D55E00") +
  labs(
    y = "MV",
    title = "Température moyenne journalière"
  )

components(dcmp_stl_tavg) %>% autoplot()

3.4.2.2 Data daily

## Desaisonnalisation de la conso corrigé avec STL 
dcmp_conso <- data_est_day %>%
  model(stl = STL(conso_corrige_dju))
components(dcmp_conso) %>%
  as_tsibble() %>%
  autoplot(conso_corrige_dju, color="gray") +
  geom_line(aes(y=trend), color = "#D55E00") +
  labs(
    y = "MV",
    title = "electricity consommation"
  )

components(dcmp_conso) %>% autoplot()

## Ajout colonne adjusted au df
data_est_day$conso_corrige_adjusted <- components(dcmp_conso)$season_adjust
## Desaisonnalisation de TAVG
dcmp_tavg <- data_est_day %>%
  model(stl = STL(TAVG))
components(dcmp_tavg) %>%
  as_tsibble() %>%
  autoplot(TAVG, color="gray") +
  geom_line(aes(y=trend+remainder), color = "#D55E00") +
  labs(
    y = "MV",
    title = "Température moyenne journalière"
  )

components(dcmp_tavg) %>% autoplot()

3.4.3 Dessaisonalisation X13

Une méthode plus récente proposé par la Spain Bank est ARIMA SEATS X13, elle s’applique uniquement aux données mensuelles

dcmp_seats_conso <- data_est_month %>%
  model(seats = X_13ARIMA_SEATS(conso_corrige_dju ~ seats())) %>%
  components()
autoplot(dcmp_seats_conso) +
  labs(title =
    "Decomposition of conso_corrigé using SEATS")

dcmp_seats_conso %>%
  ggplot(aes(x = `yearmonth(Date_only)`)) +
  geom_line(aes(y = conso_corrige_dju, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "MW",
       title = "conso_corrigé (avg/month)") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

## Ajout colonne adjusted au df
data_est_month$conso_corrige_adjusted <- dcmp_seats_conso$season_adjust
dcmp_seats_tavg <- data_est_month %>%
  model(seats = X_13ARIMA_SEATS(TAVG ~ seats())) %>%
  components()
autoplot(dcmp_seats_tavg) +
  labs(title =
    "Decomposition of Temperature (avg/month) using SEATS")

dcmp_seats_tavg %>%
  ggplot(aes(x = `yearmonth(Date_only)`)) +
  geom_line(aes(y = TAVG, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "°C ",
       title = "Temperature (avg/month)") +
  scale_colour_manual(
    values = c("gray", "#0072B2", "#D55E00"),
    breaks = c("Data", "Seasonally Adjusted", "Trend")
  )

Nous percevons bien ici l’impact du réchauffement climatique, avec une montée de 2° en moyenne sur 7 ans

3.4.4 Transformation de BoxCox

Nous pouvons utiliser une transformation de boxCox sur notre variable réponse afin normaliser les résidus. Nous choississons le parametre \(\\lambda\) grâce a l’indice de Guerrero
Cette transformation est utilisable directement dans fable, qui s’occupera alors de la transformation retour

## Transformation de Box Cox avec indice de Guerrero
lambda_guerrero_day <- data_est_month %>%
  features(Consommation, features = guerrero) %>%
  pull(lambda_guerrero)
data_est_month %>%
  autoplot(box_cox(Consommation, lambda_guerrero_day)) +
  #geom_line(aes(y=Consommation), colour = "#D55E00") +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed electricity consommation with $\\lambda$ = ",
         round(lambda_guerrero_day,2))))

gamma_guerrero_day <- data_est_month %>%
  features(TAVG, features = guerrero) %>%
  pull(lambda_guerrero)
data_est_month %>%
  autoplot(box_cox(TAVG, gamma_guerrero_day)) +
  geom_line(aes(y=TAVG), colour = "#D55E00") +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed température with $\\lambda$ = ",
         round(gamma_guerrero_day,2))))

4 Modélisation

4.1 Modèle lissage exponentielle

4.1.1 Modèle simple

# Estimate parameters
fit <- data_est_month %>%
  model(ETS(Consommation ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
  forecast(h ="1 years")
fc %>%
  autoplot(data_est_month) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="MW", title="Consommation") +
  guides(colour = FALSE)

4.1.2 Modèle double

# Estimate parameters
fit <- data_est_month %>%
  model(ETS(Consommation ~ error("A") + trend("A") + season("N")))
fc <- fit %>%
  forecast(h ="1 years")
fc %>%
  autoplot(data_est_month) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="MW", title="Consommation") +
  guides(colour = FALSE)

4.1.3 Holt-Winters

# Estimate parameters
fit <- data_est_month %>%
  model(ETS(Consommation ~ error("A") + trend("A") + season("A")))
fc <- fit %>%
  forecast(h ="1 years")
fc %>%
  autoplot(data_est_month) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="MW", title="Consommation") +
  guides(colour = FALSE)

4.1.4 Auto

Testons enfin un dernier modèle en full auto

# Estimate parameters
fit <- data_est_month %>%
  model(ETS(Consommation))
fc <- fit %>%
  forecast(h ="1 years")
report(fit)
## Series: Consommation 
## Model: ETS(M,N,M) 
##   Smoothing parameters:
##     alpha = 0.003656971 
##     gamma = 0.0002115003 
## 
##   Initial states:
##         l      s1       s2        s3        s4        s5        s6        s7
##  160457.7 1.17136 1.098812 0.9908194 0.8667747 0.7817055 0.8653705 0.8523303
##         s8        s9     s10     s11      s12
##  0.9021363 0.9519064 1.12806 1.13257 1.258155
## 
##   sigma^2:  0.0013
## 
##      AIC     AICc      BIC 
## 1839.675 1846.734 1876.137
fc %>%
  autoplot(data_est_month) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="MW", title="Consommation") +
  guides(colour = FALSE)

4.2 Modele ARIMA

AR = Auto regressive

I = Integrated

MA = Mooving average

S = Seaseonal

Approche plutôt basé sur l’autocorrélation

Nous n’utilisons plus les Trend et saisonnalité pour modéliser le processus, nous devons donc stationnarisé la série temporelle, la rendre independante du temps.

4.2.1 Stationnarisation

Nous regardons dans un premier temps si nous avons besoin de différencier grâce a un test de racine unitaire

data_est_month %>% 
  features(Consommation,unitroot_nsdiffs)
## # A tibble: 1 x 1
##   nsdiffs
##     <int>
## 1       1

Le test saisonnier nous préconise de différencier 1 fois la variable consommation en saisonnalité

Faisons le même test mais avec le test en tendance

data_est_month %>% 
features(Consommation, unitroot_ndiffs)
## # A tibble: 1 x 1
##   ndiffs
##    <int>
## 1      0

Analysons le graphes des résidus de la consommation différencié

data_est_month %>%
  gg_tsdisplay(difference((Consommation),12), plot_type='partial')

Recommençons pour la variable conso désaisonnalisé.

data_est_month %>% 
  features(conso_corrige_adjusted,unitroot_nsdiffs)
## # A tibble: 1 x 1
##   nsdiffs
##     <int>
## 1       0
data_est_month %>% 
  features(conso_corrige_adjusted,unitroot_ndiffs)
## # A tibble: 1 x 1
##   ndiffs
##    <int>
## 1      0

Le test confirme que cette variables est déjà differencié en tendance et saisonnalité

data_est_month %>% 
  gg_tsdisplay((conso_corrige_adjusted), plot_type='partial')

4.2.2 ARIMA

Nous utilisons la variable Consommation corrigé de l’effet température et désaisonnalisé

fit_arima <- data_est_month %>%
  model(
    auto = ARIMA(conso_corrige_adjusted),
    arima012011 = ARIMA(conso_corrige_adjusted ~ pdq(3, 0, 2) + PDQ(0, 0, 0))## PDQ a 0 pour indiquer qu'il n'y a pas de tendance saisonière
   # model_search = ARIMA(conso_corrige_dju ~ pdq(p=0:2, d=1, q=0:2) + PDQ(0,1,1))
    )
forecast(fit_arima,h=12) %>% 
autoplot(data_est_month)

fit_arima %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key:     Model name [2]
##   `Model name`                 Orders
##   <chr>                       <model>
## 1 auto         <ARIMA(1,0,1) w/ mean>
## 2 arima012011  <ARIMA(3,0,2) w/ mean>
glance(fit_arima) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 2 x 6
##   .model        sigma2 log_lik   AIC  AICc   BIC
##   <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 auto        6303449.   -775. 1559. 1559. 1568.
## 2 arima012011 6407252.   -774. 1563. 1564. 1580.
fit_arima %>% select(auto) %>%  report()
## Series: conso_corrige_adjusted 
## Model: ARIMA(1,0,1) w/ mean 
## 
## Coefficients:
##          ar1      ma1    constant
##       0.7702  -0.5465  29134.9501
## s.e.  0.1291   0.1586    119.0057
## 
## sigma^2 estimated as 6303449:  log likelihood=-775.32
## AIC=1558.64   AICc=1559.15   BIC=1568.37

On remarque que le modèle ne performe pas très bien du fait que notre variable ressemble de plus en plus a du bruit blanc et donc très difficile de prédire

4.2.3 SARIMA avec correction température

p=p= order of the autoregressive part;
d=d= degree of first differencing involved;
q=q= order of the moving average part. 
fit_sarima <- data_est_month %>%
  model(
    auto = ARIMA(conso_corrige_dju),
    arima012011 = ARIMA(conso_corrige_dju ~ pdq(0, 0, 2) + PDQ(0, 1, 1))
   # model_search = ARIMA(conso_corrige_dju ~ pdq(p=0:2, d=1, q=0:2) + PDQ(0,1,1))
    )
forecast(fit_sarima,h=24) %>% 
autoplot(data_est_month)

fit_sarima %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key:     Model name [2]
##   `Model name`                    Orders
##   <chr>                          <model>
## 1 auto         <ARIMA(1,0,1)(2,1,0)[12]>
## 2 arima012011               <NULL model>
glance(fit_sarima) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 1 x 6
##   .model    sigma2 log_lik   AIC  AICc   BIC
##   <chr>      <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 auto   12303309.   -691. 1391. 1392. 1402.

Nous devons maintenant réintroduire la correction de température

Nous pouvons maintenant comparé les modèle grace a la cross-validation

data_est_month_tr <- data_est_month %>%
stretch_tsibble(.init = 3, .step = 1)

fc <- data_est_month_tr %>%
    model(
    auto = ARIMA(conso_corrige_dju)

    ) %>%
  forecast(h = 24) %>%
  group_by(.id) %>%
  mutate(h = row_number()) %>%
  ungroup()
fc %>%
  accuracy(data_est_month, by = c("h", ".model")) %>%
  ggplot(aes(x = h, y = MAPE)) +
  geom_point()

4.2.4 SARIMAX (variable exogène)

Ici, nous ajoutons directement des variables éxogène à notre modèle

training <- data_est_month %>% filter(year(`yearmonth(Date_only)`) <= 2018)
test <- data_est_month %>% filter(year(`yearmonth(Date_only)`) > 2018)
fit_sarimax <- data_est_month %>%
 model(
   auto = ARIMA(Consommation~ TAVG + I(TAVG^2)),    ## Auto
   sarimax101210 = ARIMA(Consommation~ TAVG + I(TAVG^2)+pdq(1,0,1)+PDQ(2,1,0)) ## idem modele correction
   
 ) 
fit_sarimax %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key:     Model name [2]
##   `Model name`                                  Orders
##   <chr>                                        <model>
## 1 auto          <LM w/ ARIMA(0,0,1)(2,1,0)[12] errors>
## 2 sarimax101210 <LM w/ ARIMA(1,0,1)(2,1,0)[12] errors>
glance(fit_sarimax) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 2 x 6
##   .model           sigma2 log_lik   AIC  AICc   BIC
##   <chr>             <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 auto          10431952.   -685. 1382. 1383. 1396.
## 2 sarimax101210 10538566.   -685. 1383. 1385. 1399.
fit_sarimax %>% 
  select(auto) %>% 
  report()
## Series: Consommation 
## Model: LM w/ ARIMA(0,0,1)(2,1,0)[12] errors 
## 
## Coefficients:
##          ma1     sar1     sar2        TAVG  I(TAVG^2)
##       0.4284  -0.6336  -0.4129  -4795.0171   119.6446
## s.e.  0.1292   0.1225   0.1234    317.9231    14.0748
## 
## sigma^2 estimated as 10431952:  log likelihood=-685.01
## AIC=1382.03   AICc=1383.32   BIC=1395.69
fit_sarimax %>% accuracy()
## # A tibble: 2 x 10
##   .model        .type       ME  RMSE   MAE    MPE  MAPE  MASE RMSSE     ACF1
##   <chr>         <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>
## 1 auto          Training -360. 2885. 2087. -0.254  1.31 0.297 0.286 -0.0418 
## 2 sarimax101210 Training -386. 2878. 2105. -0.272  1.33 0.300 0.285 -0.00312
# Time series cross-validation accuracy
data_est_month_tr <- data_est_month %>%
  stretch_tsibble(.init = 1, .step = 12)
data_est_month_tr
## # A tsibble: 259 x 14 [1M]
## # Key:       .id [7]
##    `yearmonth(Date_only)`   TAVG  TMAX   TMIN Consommation Renouvelable
##                     <mth>  <dbl> <dbl>  <dbl>        <dbl>        <dbl>
##  1             2013 janv.  2.25   4.33  0.103      202862.       42650.
##  2             2013 janv.  2.25   4.33  0.103      202862.       42650.
##  3             2013 févr.  0.968  3.71 -1.45       191247.       41756.
##  4              2013 mars  3.56   7.56  0.103      195760.       40565.
##  5              2013 avr. 10.9   16.0   6.30       161972.       48322.
##  6               2013 mai 12.9   17.2   9.05       150792.       48971.
##  7              2013 juin 18.7   24.0  13.2        137521.       48742.
##  8             2013 juil. 22.4   28.3  15.9        140700.       46800.
##  9              2013 août 19.7   26.2  13.8        124791.       34938 
## 10             2013 sept. 15.8   21.5  11.3        141255.       37687.
## # ... with 249 more rows, and 8 more variables: dju_calc_chauffage <dbl>,
## #   dju_calc_clim <dbl>, dju_calc_total <dbl>, PRCP <dbl>,
## #   conso_corrige_dju <dbl>, sum(conso_corrige_temp) <dbl>,
## #   conso_corrige_adjusted <dbl>, .id <int>
# TSCV accuracy
fit <- data_est_month_tr %>%
  model(ARIMA(Consommation)) %>%
  forecast(h = "2 years") %>%
  accuracy(data_est_month)
forecast(fit_sarimax,new_data = test) %>%
  autoplot(data_est_month) +
  labs(title = "Electricity demand",
       y = "MW") + coord_cartesian(xlim = c(as_date("2018-01-01"),as_date("2019-12-31")))

4.3 Modèle deep learning

Cette méthode est basé sur les réseaux de neurones artificiels, nous pouvons aussi y ajouter des variables explicatives.

Le temps de calcul étant relativement long, nous travaillons sur le data month. Vu que nous utilisons une variable exogene dans notre modèle, nous splittons le dataset en Train/Test

fit_mnn <- training %>% 
  model(
        AutoX = NNETAR(Consommation , xreg = training$TAVG),
        Auto = NNETAR(Consommation)
        )
forecast(fit_mnn,new_data = test, times = 100) %>% ## times = 200 permets de limiter le temps de calcul en dimuniant la précisions des intervalles de confiance
  autoplot(data_est_month) +
  labs(title = "Electricity demand",
  y = "MW") + coord_cartesian(xlim = c(as_date("2018-01-01"),as_date("2019-12-31")))  ## Restiction pour ne voir que 2 ans de 2018 a 2019

accuracy(fit_mnn)
## # A tibble: 2 x 10
##   .model .type       ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AutoX  Training -4.62 7281. 5219. -0.191  3.09 0.747 0.715 0.348
## 2 Auto   Training  8.03 7366. 5237. -0.186  3.10 0.750 0.724 0.334

Les résultats sont convainquant avec un MAPE de 3%

4.4 Modele Prophet

Le modèle Prohet est un modèle récent de 2018 (Taylor and Letham 2018) viens directement de Facebook, basé sur des régressions avec des séries de Fourier.

Ses avantage sont sa simplicité d’utilisation et son intégration des effets vacances, l’incovenient principale est que nous ne maitrisons pas les choix de modèle, basé sur des critère Bayésien

library(prophet)
## Loading required package: Rcpp
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
## 
##     %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
##     flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
##     splice
library(dygraphs)

Nous importons un fichiers avec les jours fériés Francais, nous pourrions importer les dates des vacances mais cela nécessite plus de manipulation. De plus, ils devrais être integrés prochainement directement dans Prophet

jours_feries <- read.csv("data/holidays_fra.csv",fileEncoding="UTF-8")

Nous devons renomer les variables date (ds) et d’interet (y)

data_prophet_day <- data_est_day %>%
  rename(ds = Date_only, y = Consommation) 
data_prophet_month <- data_est_month %>%
  rename(ds = `yearmonth(Date_only)`, y = Consommation)

Nous renseignons le modele

m<- prophet(holidays = jours_feries) ## ajout des jours feries
m<- add_regressor(m,'TAVG') ## ajout variables expliquatives
m<- fit.prophet(m,data_prophet_day) ## fit du model
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Nous devons crées un dataframme future a partir duquel nous pourrons faire les prévision. Pour le moment, nous avons que les dates dedans, il nous faut les prévisions de la variable exliquative. Nous pourrions l’importer depuis un organisme météo, mais j’ai choisi de la prévoir directement à l’aide de Prophet

data_prophet_temp2 <- data_prophet_day %>%
  rename(Consomation=y,y = TAVG) ## renommage des variables
m2 <- prophet(data_prophet_temp2) ## on renseigne le df 
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
future <- make_future_dataframe(m2, periods = 365) ## on crée un df future de 1 an
forecast <- predict(m2, future) ## on realise les predictions
forecast_temp_2020 <- data.frame(forecast[c('ds', 'yhat')]) %>%  ## on transforme les prédictions en df
  rename(TAVG=yhat) ## on renome TAVG

Nous pouvons maintenant lancer les prévisions.

forecast <- predict(m, forecast_temp_2020) ## predictions

prophet_plot_components(m, forecast)## composante

dyplot.prophet(m, forecast,uncertainty = TRUE)## Affichage dynamique

Vérifions maintenant la qualité des prédictions à l’aide d’une cross validation

df.cv <- cross_validation(m, horizon=365, units='days')
## Making 7 forecasts with cutoffs between 2016-01-01 and 2018-12-31
plot_cross_validation_metric(df.cv, metric='mape')

Nous obtenons un résultat très intéressant, avec une erreur moyenne de 5%

Conclusion

Taylor, Sean J., and Benjamin Letham. 2018. “Forecasting at Scale.” The American Statistician 72 (1): 37–45. https://doi.org/10.1080/00031305.2017.1380080.